Specifically, we will explore the results of 54 datapoints from 4 studies conducted in 1 ocean on 2 species within 2 genera.
## n
## 1 54
## Ref n
## 1 DS19 38
## 2 DS01 8
## 3 DS06 6
## 4 DS02 2
## Ref_name n
## 1 Ricardo et al. (2017) 38
## 2 Babcock and Davies (1991) 8
## 3 Hodgson (1990a) 6
## 4 Babcock and Smith (2000) 2
## Ocean n
## 1 Pacific 54
## Gsp n
## 1 Acropora_millepora 48
## 2 Pocillopora_damicornis 6
## Updated_Genus n
## 1 Acropora 48
## 2 Pocillopora 6
For larval settlement, we will explore the results of 30 datapoints from 7 studies conducted in 2 oceans on 4 species within 3 genera.
## n
## 1 30
## Ref n
## 1 SS13b 5
## 2 SS13c 5
## 3 SS13d 5
## 4 SS23a 4
## 5 SS24a 4
## 6 SS24b 4
## 7 SS11c 3
## Ref_name n
## 1 Humanes et al. 2017b 15
## 2 Te 1992 8
## 3 Rushmore 2016 Chapter 2 4
## 4 Gilmour 1999 3
## Ocean n
## 1 Pacific 26
## 2 Atlantic 4
## Gsp n
## 1 Acropora_tenuis 15
## 2 Pocillopora_damicornis 8
## 3 Porites_astreoides 4
## 4 Acropora_digitifera 3
## Updated_Genus n
## 1 Acropora 18
## 2 Pocillopora 8
## 3 Porites 4
For larval survival, we will explore the results of 63 datapoints from 4 studies conducted in 1 ocean on 5 species within 2 genera.
## n
## 1 63
## Ref n
## 1 SS20d 37
## 2 SS13b 5
## 3 SS13c 5
## 4 SS13d 5
## 5 SS24a 4
## 6 SS24b 4
## 7 SS11b 3
## Ref_name n
## 1 Ricardo et al. 2016 37
## 2 Humanes et al. 2017b 15
## 3 Te 1992 8
## 4 Gilmour 1999 3
## Ocean n
## 1 Pacific 63
## Gsp n
## 1 Acropora_tenuis 31
## 2 Acropora_millepora 16
## 3 Pocillopora_damicornis 8
## 4 Pocillopora_acuta 5
## 5 Acropora_digitifera 3
## Updated_Genus n
## 1 Acropora 50
## 2 Pocillopora 13
First let’s explore all data from all species for which there is data about ‘limited settlement’ and ‘reduced larval survival’ as a result of exposure to sediment.
Now let’s calculate the thresholds, based on the binary data explored above.
## LOAEL
## 1 1
## NOAEL
## 1 1
## LOAEL
## 1 57.79
## NOAEL
## 1 34.6
## LOAEL
## 1 30
## NOAEL
## 1 29.5
## n
## 1 45
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: settleDS2
## Models:
## modDS4: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp) +
## modDS10: (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS4 3 41.580 47.000 -17.790 35.580
## modDS7 3 41.515 46.935 -17.758 35.515 0.0651 0 <2e-16 ***
## modDS10 4 43.515 50.742 -17.758 35.515 0.0000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## df AIC
## modDS1 2 39.58049
## modDS7 3 41.51537
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## Data: settleDS2
##
## AIC BIC logLik deviance df.resid
## 41.5 46.9 -17.8 35.5 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7321 -0.7847 0.0016 0.4885 1.2405
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 1.089 1.044
## Number of obs: 45, groups: Ref, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.06776 0.94288 -0.072 0.9427
## Sed_level_stand_mg 0.07684 0.04664 1.648 0.0994 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 0.8444444
##
## p 0 1
## 0 9 4
## 1 3 29
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9167
## R2m R2c
## theoretical 0.9847224 0.9885223
## delta 0.9765060 0.9802742
## $Ref
## Groups Name Std.Dev.
## Ref (Intercept) 1.0436
## [1] 2.839528
## Est LL UL
## (Intercept) 0.9541154 0.2650232 3.434930
## Sed_level_stand_mg 1.0547038 0.9899510 1.123692
There is suggestive, but non-significant, evidence that for every doubling of exposure concentration of deposited sediment, the odds of reduced settlement rate increase by 1.05 times (95% CI 0.99, 1.12; GLMM z = 1.648, p = 0.09) after accounting for variability among studies.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS7 <- predict(modDS7, newdata=settleDS2, type="response")
settleDS3 <- cbind(settleDS2, pred_modDS7)
settleDS_plot <- ggplot(data = settleDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_limited_settlement,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Reduced settlement rate due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(0.1,1100),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS7), inherit.aes=FALSE) +
theme_bw()
settleDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(settleDS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
settleDS2$Sed_level_stand_mg <- j
predict(modDS7, newdata = settleDS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
settleDS_plot3 <- ggplot() +
geom_point(data = settleDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_limited_settlement,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of reduced settlement\nrate due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,1100), breaks=c(0.01,0.1,1,10,100,1000,loaelDS),
label=c("0.01","0.1","1","10","100","1000",round(loaelDS, digits=1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelDS, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
settleDS_plot3
## n
## 1 20
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: settleSS2
## Models:
## modSS4: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp)
## modSS7: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## modSS10: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Gsp) +
## modSS10: (1 | Ref)
## modSS13: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 20.127 23.114 -7.0636 14.127
## modSS7 3 19.416 22.403 -6.7079 13.416 0.7115 0 <2e-16 ***
## modSS10 4 22.127 26.110 -7.0636 14.127 0.0000 1 1
## modSS13 4 22.127 26.110 -7.0636 14.127 0.0000 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## df AIC
## modSS1 2 20.75610
## modSS7 3 19.41575
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_limited_settlement ~ Sed_level_stand_mg + (1 | Ref)
## Data: settleSS2
##
## AIC BIC logLik deviance df.resid
## 19.4 22.4 -6.7 13.4 17
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.55090 -0.01777 -0.01380 -0.01347 1.72158
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 120.9 10.99
## Number of obs: 20, groups: Ref, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.522305 5.925743 -1.438 0.150
## Sed_level_stand_mg 0.001455 0.005675 0.256 0.798
## convergence code: 0
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_limited_settlement ~ log10_conc + (1 | Ref)
## Data: settleSS2
##
## AIC BIC logLik deviance df.resid
## 14.3 17.3 -4.1 8.3 17
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.01213 0.00000 0.00000 0.00000 0.05431
##
## Random effects:
## Groups Name Variance Std.Dev.
## Ref (Intercept) 18037 134.3
## Number of obs: 20, groups: Ref, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -95.32 46.01 -2.072 0.0383 *
## log10_conc 27.27 14.22 1.918 0.0551 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 1
##
## p 0 1
## 0 17 0
## 1 0 3
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 1
## R2m R2c
## theoretical 1.781484e-02 9.998209e-01
## delta 7.265477e-14 4.077598e-12
## $Ref
## Groups Name Std.Dev.
## Ref (Intercept) 134.3
## [1] 2.118497e+58
## Est LL UL
## (Intercept) 4.001764e-42 2.730730e-81 5.864407e-03
## log10_conc 6.978395e+11 5.493766e-01 8.864228e+23
There is suggestive, though non-significant, evidence that for every 10-fold increase in exposure concentration of suspended sediment, the odds of reduced settlement rate increase (GLMM z = 1.918, p = 0.055), while accounting for variability among studies.
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS7_log <- predict(modSS7_log, newdata=settleSS2, type="response")
settleSS3 <- cbind(settleSS2, pred_modSS7_log)
settleSS_plot <- ggplot(data = settleSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_limited_settlement,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Reduced settlement rate\ndue to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS7_log), inherit.aes=FALSE) +
theme_bw()
settleSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(settleSS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
settleSS2$log10_conc <- j
predict(modSS7_log, newdata = settleSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
plotdat <- plotdat %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_mg_L", "mg_L")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = log10_mg_L)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = log10_mg_L)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = mg_L, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
settleSS_plot2 <- ggplot() +
geom_point(data = settleSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_limited_settlement,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of reduced settlement\nrate due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,max(settleSS2$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS1),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS1,digits = 1))) +
geom_ribbon(data = plotdat, aes(x = mg_L, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = mg_L, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = mg_L, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS1, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
settleSS_plot2
## n
## 1 52
## df AIC
## modSSb1 2 44.75628
## modSSb1b 3 40.72043
## Data: survivalSS2
## Models:
## modSSb4: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Gsp)
## modSSb7: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Ref)
## modSSb10: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSSb10: Ref)
## modSSb13: Binary_larval_survival ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSSb4 3 44.117 49.971 -19.059 38.117
## modSSb7 3 43.904 49.758 -18.952 37.904 0.2132 0 <2e-16 ***
## modSSb10 4 44.063 51.868 -18.031 36.063 1.8415 1 0.1748
## modSSb13 4 45.695 53.500 -18.848 37.695 0.0000 0 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## df AIC
## modSSb1 2 44.75628
## modSSb4 3 44.11739
##
## Call:
## glm(formula = Binary_larval_survival ~ Sed_level_stand_mg, family = binomial(link = "logit"),
## data = survivalSS2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.7033 -0.5272 -0.5016 -0.4975 2.0679
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.0356448 0.5215749 -3.903 9.51e-05 ***
## Sed_level_stand_mg 0.0007646 0.0012874 0.594 0.553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 41.087 on 51 degrees of freedom
## Residual deviance: 40.756 on 50 degrees of freedom
## AIC: 44.756
##
## Number of Fisher Scoring iterations: 4
## [1] 0.8653846
##
## p 0 1
## 0 45 7
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.6429
## R2m R2c
## theoretical 0.014274072 0.014274072
## delta 0.005519128 0.005519128
There is no significant relationship between exposure concentration of suspended sediment and the odds of larval mortality (GLMM z = 0.762, p = 0.446).
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb1 <- predict(modSSb1, newdata=survivalSS2, type="response")
survivalSS3 <- cbind(survivalSS2, pred_modSSb1)
survivalSS_plot <- ggplot(data = survivalSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_larval_survival,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Larval mortality due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(1,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb1), inherit.aes=FALSE) +
theme_bw()
survivalSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(survivalSS2, seq(from = min(Sed_level_stand_mg), to = max(Sed_level_stand_mg), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
survivalSS2$Sed_level_stand_mg <- j
predict(modSSb1, newdata = survivalSS2, type = "response")
})
# get the means with lower and upper quartiles
plotdat <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_level_stand_mg values and convert to data frame
plotdat <- as.data.frame(cbind(plotdat, jvalues))
# better names and show the first few rows
colnames(plotdat) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_conc")
#head(plotdat)
# plot average marginal predicted probabilities
ggplot(plotdat, aes(x = Sed_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat, aes(x = Sed_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
survivalSS_plot2 <- ggplot() +
geom_point(data = survivalSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_larval_survival,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of larval\nmortality due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(1,max(survivalSS2$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits = 1))) +
geom_ribbon(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
survivalSS_plot2